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LIST OF SYMBOLS AND ABBREVIATIONS 



A Arbitrary constant in the stream function 
a Earth's radius 

F Zonal flux term 

g U Flux term in ^-equation of motion 

g V Flux term in ^-equation of motion 

G Meridional flux term 

I Number of points around latitude circle 

i Grid index in the ^-direction 

j Grid index in the ^-direction 

k Vertical grid index 

m 1/ (a cos 4>) 

mb Millibars 

NACA National Advisory Committee on Aeronautics 
n 1/a 

P Pole on the ij index system 

S Area-pressure weighted vertical velocity 

t Time 

u Zonal wind 

v Meridional wind 

n Meridional coordinate of the curvilinear coordinate system 

Ap Distance increment in the meridional direction 
X Longitude 

^ Zonal coordinate of the curvilinear coordinate system 
Distance increment in zonal direction 
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rr Terrain pressure 
\\ Area-weighted terrain pressure 

a Dimensionless vertical coordinate 

Aa Vertical increment in the sigma coordinate system 
a Measure of vertical velocity 
Latitude 

0 Angulaf velocity of the earth 
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I. INTRODUCTION 



Most global primitive equation models now in use employ spherical 
coordinates (see Haltiner and Williams [ 1975]). Flow .crossing the pole 
presents a potential problem in this coordinate system. Mihok and Kaitala 
(1976) have described a global prediction model which is in the last stages 
of development at the Fleet Numerical Prediction Central. McCollough (1974) 
in testing an earlier version of this model observed the development of large 
gradients in the surface pressure near the pole when a particular real data 
set was used. Maher (1974) examined this problem more closely by using an- 
alytic initial data which gave a strong flow across the pole. His solutions 
showed that the pressure field was spuriously disturbed by the finite dif- 
ferencing near the pole. He was able to reduce the effect by various types 
of smoothing. 

In this report we will test a somewhat different finite difference 
scheme and we will test a procedure for controlling the problem. The Navy 
Environmental Prediction Research Facility is now testing the global pre- 
diction model which was developed by Arakawa and Mintz (1974). In this 
study we will employ the model which was described by Monaco and Williams 
(1975). This model is a slightly simplified version of the Arakawa model 
and it differs from the FNWC model mainly with respect to the spatial 
staggering of the variables. 
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II. FINITE DIFFERENCE EQUATIONS 



The basic differential and difference equations used ip. this study 
are given in Monaco and Williams (1975).^ Here we will reproduce only 
those difference equations which have been changed. Fig. 1 shows the 
arrangement of variables in the model. The list of symbols gives the 
various definitions. The finite difference continuity equation has 
the following form: 





0 . 



( 2 . 1 ) 



The mass fluxes are defined 



k 




( 2 . 2 ) 




(2.3) 



where 





The bar in (2.2) represents zonal smoothing as described inMW 
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Hereafter we will refer to this reference as MW 
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! 





RED 

ION 



Figure 1. An example of the notation used to describe the finite dif- 
ference equations. The continuity equation is described on a "ir-centered" 
grid in which the F and G symbols are flux calculations in their re- 
spective directions. The "u-centered" grid is an example used to de- 
scribe the ^-component of the equation of motion where F U and g J are 
also flux calculations. 
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u-CENTE 

NOTAT 



The left hand side of the zonal momentum equation is: 



+4tF i+i,j <u i+i,j + u ij : 



„u . , N k , u , , N ku / . k, 

- F. A a(u. 4 + u ,- , a> + §,• i 4 Jk( u n - a 4.1 + u ,- a) - g-: A_i( u i a + u ,- a J 1 

1- "2 9 J 1 I 5 J 1 ? J 1 j J^l 2 1 * J 1 



1 i r *u,k-R , k+2 , k x *u,k-l f k k-2xi 

+ — r iK S. ’ . (u. . + u # . S. 5 . (u. . + u. .) J 

Ac k i,J i>J i,J i,J i,J 



(2.4) 



where 



*' i, j 8^i+|,j+l + lli-|,j+l + + 1|i+|,j-l +2 ^li+§,j + l\i-§,j^’ (2 ’ 5) 



•u i r ; 



S7 . = ^[S., ... + S. 1 ... + S. 1 . . + S . - . . + 2(3., . + S. a. .)]. (2.6) 

i,j 8 1+2, j +1 1-2, J +1 1-2, J-l i+ 5 ,J“l i-2»J 



The above definitions (2.5) and ( 2 . 6 ) are more complicated than those in MW 
and they are the same as those used by Arakawa and Mintz (1974). The quan- 
tities F U and g U are defined as follows: 



F U ^i . s 7 (F*i ... + 2F*i . + F* , . ) , 

1+ 2jJ 4 1+ 5>J + 1 l+g+J 



(2.7) 






= 7~ (G ■ * . l • + Gy ,i . . + G^ i . + Gt i . .) , 
4 1 + 2 , j i+f,J+l i-2,J i-5, J+1 



( 2 . 8 ) 



where 



F* . = i(F. , . + F. a. .) , 

1+2 » J 1-5, J 



(2.9) 



G* . — g’ (G . . 1 + G. . i) 

i,J i,J+5 i,J-5 



( 2 . 10 ) 
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The left hand side of the meridional momentum equation takes the 
same form as (2,4) when u is replaced by v. The other quantities are 
defined as follows: 

(2.13) 

(2.14) 

The definitions (2.13) and (2.14) are more complicated than those in MW 
and they are the same as those used by Arakawa. 

Arakawa and Mintz (1974) derived the expressions (2.5) and (2.6) for 
IT and S at the u points in the following way. The u field in (2.4) is set 
equal to a constant. Then the continuity equation (2.1) is averaged from 
surrounding points in such a way that the average has the same form as 
(2.4) with u = constant. This is a reasonable requirement and it is nec- 
essary for energy conservation. This average gives the expressions (2.5) 
and (2.6). The definitions (2.13) and (2.14) are derived in the same way 
by setting v = constant in the appropriate equations. 

These difference expressions must be modified near the pole because 
some of the quantities are not defined at the pole. Fig. 2 shows the 
points near the pole which are used to evaluate the pressure change at 
the pole. The pole is composed of I points which are denoted by (i, P) . 



F V i = — (F* *1 + F* i 4- F* i 4- F* i . 

i+i, j i+1 s j+i j+l i,j-§ i+l, j“t 



( 2 . 11 ) 



S i» j+2 ^ ^ G i+1, j+i + 2G i , j+2 + G i“l, j+2^ ’ 



( 2 . 12 ) 



V i, j 8^Ti+i , j+i + jli-l, j+i + /ii-l, j~i + ||i+i, j~i + 2< ^i, j+i + |U, j-|^ 



§ 1 * • • • ** 

s v . = S . , , . , 1 +S . , + S. , . ! +S . , , . 1 + 2 (S . .^1 + S. . i) ] 

i,j 8 1 + 1 , j+t i-l, J+i i-l, j-i i+l, j-t i , J+2 i,J-2 
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n 




i 



Figure 2. Each index i,P in the continuity equation is 

represented by the shaded area, tt at the poles 
can only change as a result of G. The thermo- 
dynamic equation and vertical velocity are treated 
in the same manner. 
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The continuity equation applied at each of these points is written: 




(2.15) 



After each time step the surface pressure at the pole (P) is set equal 
to the average: 



If we sum (2.18) in the vertical and use (2.16) we can see that the pres- 
sure at the pole changes in proportion to the net mass flux across the 
latitude circle at j = P-§. 

Fig. 3 shows the quantities which are needed to predict u at j = P-1. 
The left hand side of the zonal momentum equation at j = P-1, takes the form 




(2.16) 





+ u 





(2.17) 



where 




P-1 




(2.18) 
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*> 7 r 




Figure 3. The polar modification of the u equation of motion- 
Z 7 and g are flux terms. The shaded portion 
represents the area associated with each variable 
u . 
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K 



~ .3, 



P-1 



ilp + + /]n J k,P-l ) + sffi-i. 



P-2 



+ 



II 



i-ti,P-2 } 



(2.19) 



S u 



i,P-l 



= 



3 * 

+ 4(s 



8 V i-#,P-l 



+ S 



i+2 ,P-1 



) + 4(S. 



^ P-2 

, , r z. 



+ S. 



|,P-2 ) * 



(2.20) 



The other quantities have the same definitions as above. The relations 
(2.19) and (2.20) are different from those used by MW and they correspond 
to the relations derived by Arakawa. 

Fig. 4 shows the quantities which are needed to predict v at j = P-g-. 

The left hand side of the meridional momentum equation at j = P-^ 3 takes the form 



dt (r7 v } i,p-§ + 2 [F i+|,p4 (v i+i,p-i + v i,p4 } 



~ F i-i,P-! (V i,P4 + V i-l,P4 } " 8 i,P-l (V i,P4 + V i,P-3/2 } 1 

+ 1 ^|[S V ’ k+ Vfp 1 +vj x ) - S V ’ k_ V k p i) +v k “p_i)l, (2.21) 



where 



F ?4.p-i“i (F J-i.p-i + r i.p-i )> 



( 2 . 22 ) 



7*1, P-f "^P + 8^1-1, P-1 + TTi+l,p-l + 2 Wi,P-l )5 (2.23) 



S i,P-§ “ S P + 8 (S i-l,P-l + S i+l,P-l +2S i,P-l ) * 



(2.24) 
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i 



Figure 4. The polar modification of the v equation of 

*V *V 

motion. F and g are flux terms. The shaded 
portion represents the area associated with each 
variable v. 
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The relations (2.23) and (2.24) are different from those used by MW and they 
are equal to the relations derived by Arakawa. 



The expressions for Jj 11 and Jf V near the pole were derived by Arakawa 



and Mintz (1974) in the same manner as for other latitudes. For example, 
u is set equal to a constant in Eq.(2.17) and the continuity equation is 
averaged to achieve the same form. However, the derivation used (2.15), 
but does not use the average condition (2.1$). 

III. NUMERICAL SOLUTION WITH FLOW OVER POLE 

In this section we will examine a numerical solution which was ob- 
tained with the difference equation described in section 2. The initial 
conditions, which are similar to those used by Holloway, Spelman and 
Manabe (1973), contain strong flow across the poles. The initial con- 
ditions are taken from the barotropic Rossby wave solution developed by 
Haurwitz (1940) . These initial data lead to an exact solution for baro- 
tropic horizontal motion (Neamtan [1946]). However, in our model the di- 
vergence will not remain zero, but the exact solution should be close to 
the Haurwitz solution. 

The initial zonal wind is given by 




(3.1) 



and the meridional component is given by 




(3.2) 
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These fields describe wavenumber 1 planetary flow with no mean current. 

The initial surface pressure is given by 

rr(X,*,t) = tt + K[ (■—) 2 (2cos^<t> - cos 2 * -2) 
ave 2a 

+ ^cos 4> (5 — 4 cos 2 <t>) sin A. + (~)^ cos 2 4>(3 - 2cos 2 <f>) cos \ ] . (3.3) 

3 2a 

Here TT ave is the mean surface pressure and K is a constant of proportion- 
ality. This solution for the pressure field was obtained by Phillips (1959) 
and is equivalent to the solution of the nonlinear balance equation. The 
initial temperature field is a function of height only, and it is given by 
the NACA standard atmosphere* 

Fig. 5 shows the initial surface pressure field in the vicinity of the 
pole. Note the strong geostrophic flow directly across the pole. 

The tests described in this report were carried out with a two-level 
model, with 45 points between the poles and 10 points in the east-west 
direction. A time step of 6 minutes was used. No heating, friction or 
surface topography were included. 

Fig. 6 shows the surface pressure prediction of t = 54 hours which was 
made from the above initial conditions. The basic pattern has rotated con- 
siderably from its initial orientation because the wave number 1 initial 
state excites a Rossby wave which moves rapidly westward. However, in the 
vicinity of the pole the isobars have become quite distorted. In the non- 
divergent solution obtained by Neamtan (1946) the initial disturbance rotates 
at constant angular velocity without distortion. Although in our model the 
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Figure 5. The initial surface pressure field in the vicinity of the pole 
The location of the pole is indicated by the small circle in the center. 
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Figure 6. The surface pressure field in the vicinity of the pole at 
t = 54 hours. 
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divergence is not zero, it is expected that the general behavior should 
be similar to the behavior of the nondivergent solutions. Thus the dis- 
tortions which appear in Fig. 6 near the pole are of a computational 
nature. In fact, the numerical solution ,f blew-up” at about t = 75 hours. 

In order to consider this problem further, let us examine the 
component of the velocity near the pole. In the upper part of Fig. 7 is 
the initial v component at $ = 88°; the v component at ♦ = 84° is almost 
identical. In the lower portion of Fig. 7 are the v^ fields at the 
2 latitudes for t = 54 hours. The fields are somewhat out of phase with 
each other and the component nearest to the pole is greatly amplified. 

After t = 54 hours the maximum v component continues to grow until the num- 
erical solution is completely destroyed. 



IV. CONTROL OF POLAR PROBLEM 

In Fig. 7 it was seen that the v field along the row of points next to 
the pole becomes very large and irregular as the numerical solution becomes 
unrealistic. It appears that the geostrophic adjustment process cannot 
operate effectively along the row of points next to the pole = 88° in 
our experiments). For example, if the v component is too large at one 
point, then it would modify the polar value of the pressure, which would 
change the pressure gradient at the original points. However, this pro- 
cess cannot operate freely because the polar pressure changes in propor- 
tion to all the points of <J> = 88°, not just one (see Eqs. (2.15) and (2.16)). 
It is true that various quantities near the pole such as (2.22) and (2.23) 
and (2.24) were derived by assuming consistency between the momentum change 
at <f> = 88° and the pressure changes at surrounding points. But this proof 
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A — 

2 

Figure 7. The v field as a function of longitude at t = 0 and t = 54 hours 
at the latitudes indicated 0 
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is not strictly valid because it assumes different values of tt- 

i , r 

(see Eq. (2.15)) for different longitudes whereas (2.16) is enforced at 
every time step. 

In order to avoid this problem it seems appropriate to keep only that 
portion of the v field at = 88° that can achieve geostrophic balance. This 
restricts the v field to just wave numbers 0 and 1. With wave number 0 we 
have the proper interaction with the polar pressure, since v is independent 
of longitude. The velocity vector is the same on both sides of the pole with 
wavenumber 1, so that the polar pressure should not change. Thus after every 
time step we set 



v 



i,p4 ' i i Si v i,p-i + Ksi v i,p-i cos( ¥ i) > c » s (^r> 



+ I ( Jl v i,p-i sin( ^ ) 



(4.1) 



Fig. 8 shows the surface pressure field at t = 54 hours when (4.1) is 
used every time step. In this case the isobars are still distorted near 
the pole, but not nearly as much as in Fig. 6. In Fig. 9 we see the pre- 
dicted v 1 fields at the 2 highest latitudes for t = 54 hours. These 
fields are much closer to the initial fields than those shown in Fig. 7. 
There is a small increase in the velocity at 4> = 88° and both fields are 
quite smooth. Of course v. i is forced to be smooth by (4.1). Fig. 10 
shows the surface pressure field at t = 108 hours. It is actually smoother 
than at t = 54 hours, although the later solution does have some tilt. It 
should be noted that the total disturbances amplitude has decreased by about 
20 $ at t = 108 hours. This is apparently due to the Euler -backward restart 



20 




Figure 8. The surface pressure field at t = 54 hours obtained with the 
modified v structure. 
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A-— 



2 

Figure 9. The v field as a function of longitude obtained with the 
modified v structure for t = 54 hours and t = 108 hours 0 
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Figure 10. The surface pressure field at t = 108 hours obtained with 
the modified v structure. 
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which is used in the model every 30 minutes and to the high frequency of 
the Rossby wave in this case. The v^ field near the pole is smooth and 
has a small amplitude as may be seen in the lower portion of Fig. 9. 



V. CONCLUSIONS 

It has been shown that the global model developed by Monaco and 
Williams (1975) is inaccurate when there is flow over the pole, and in 
fact the model may become unstable. This instability was attributed to 
the difficulty in achieving geostrophic adjustment near the pole. This 
is because the polar pressure is changed by the average mass flux across 
the row of points nearest to the pole. Thus the polar pressure reacts 
only weakly to the v at a single point adjacent to the pole. 

The instability was eliminated when the v. i field was replaced 

i,P - 2 

by its average plus the first wave in longitude. The resulting integra- 
tions were stable and the fields were smooth, although some distortion still 
occurred near the pole. Sadourny (1975) has formulated a global model in 
cylindrical coordinates which conserves mean square potential vorticity. 

His solutions with wave number 1 appear to be quite good, although he does 
not show the detail near the pole. 

The technique developed in this report may be useful for controlling 
the polar problem in other global models. It could be used in the model 
described by Mihok and Kaitala (1976) which is now under development at 
the Fleet Numerical Weather Central. In this application the fields u, v, 
rr and T would all have to be treated since they all occur on the same 
latitude circle. 
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